Copyright © 1973 American Telephone and Telegraph Company
The Bell System Technical Jocbnal
Vol. 52, No. 6, July-August, 1973
Printed in U.S.A.
Charge Distribution in Buried- Channel
Charge- Coupled Devices
By W. H. KENT
(Manuscript received January 11, 1973)
This paper studies charge distribution in buried-channel charge-
coupled devices. Detailed development of a one-dimensional electrostatic
model is presented and a numerical solution of the resulting nonlinear po-
tential equations is described. Graphical results show the charge-filling
mechanism and the relationship between the oxide-semiconductor interface
potential and total free positive charge.
I. INTRODUCTION
This paper describes a numerical determination of the distribution
of charge in a one-dimensional model of a buried-channel charge-
coupled device (CCD). 1 Several calculations have recently been made
of the static potential in CCD's in the absence of stored charge. 2-4 In
addition, the motion of the stored charge under dynamic conditions
has been studied by means of essentially one-dimensional models which
do not involve a true knowledge of the distribution of stored charge or
the charge-carrying capacity of the CCD. 2 - 5 However, so far it has not
been possible to calculate even the static stored-charge distribution in
a two-dimensional model of a buried-channel CCD, much less to follow
the motion of this charge under dynamic conditions. In this paper a
start is made on the problem by calculating the static distribution of
stored charge in a one-dimensional model of a buried-channel CCD.
The resulting information on the charge distribution is of interest in
itself. However, an additional important objective has been to find
numerical techniques which can be extended to the two-dimensional
problem.
The paper is divided into three parts. Section I treats the physics of
the model and Section III gives numerical results. Section II deals
briefly with the numerical techniques used and it may be omitted by
the reader without loss of continuity.
1009
< I
1010 THE BELL SYSTEM TECHNICAL JOURNAL, JULY-AUGUST 1$73
-SILICON DIOXIDE
METAL
i h -
p-TYPE SILICON
n-T'YPE SILICON
Fig. 1 — Buried-channel device structure.
The buried-channel device has a layered structure wjiich will be
modeled in one dimension as follows (Fig. 1): ' .
£ * ^ hi', silicon dioxide (Si0 2 ) with relative dielectric constant
ei/e =s 4. (e D is capacitivity of free space.)
hi ^ x ^ h 2 : p-type silicon with acceptor number density Wa(x)
and relative dielectric constant e 2 /t = 12. :
h 2 ^ x ^ oo : n-type silicon uniformly doped with constant donor
number density Nt> and dielectric constant e 3 — *i.
The point x = is a perfectly conducting boundary held at a potential
V . The potentials 4>i(x), faix), and 4>z(x), in the-Si0 2 , p-type\ and
n-type regions, respectively, satisfy the one-dimensional Poisfcon
equation,
d*<t>i{x)
= -p.oo/e.v t = 1,2,3, * ;•; \<i):
dx 2
where the pi(x), the lineal charge densities, are nonlinear functions of
the potentials <tn{x). ; ,.' ■
The functional forms of the pi(x) are determined by' the assumptions:
(i) The Si0 2 is a perfect insulator.
(ii) The generation and recombinatioh rates for hole's and electrons
are zero. > • .
(in) Hole and electron currents are zero at the time of observation.
(iv) The flat-band voltage is zero.
(v) The injected free holes and electrons are separately in
equilibrium. ' . j
Conditions (i) through (v) define a static device for which the p.-(af) are:
Pi(z) =0, "I ■ . i
P2(x) =qlp(x)-N A (xn <\: ■ (2)
Pi(*) = <£#» + *»(*) T *MlJ .
where w(x') and pCx) are the number densities df free electron^ and
holes, respectively. The most general expression' for rtt» 3B
n(z) = /"" gimrXfydti,^ ' .' , • (3)
BURIED-CHANNEL CCD 1011
where g{E) is the density of states, F e (E) is the Fermi distribution for
electrons, and E c {x) is the conduction band edge. Substitution for g(E)
and F.(E) yields
n{X) * V .. <., 1 + exp «E - E,)/kT) (4)
This specific choice for the density of states corresponds to the simplest
possible band structure. E/ is the equilibrium Fermi level for electrons
and N° c is the constant
N° c = 47r(2m c /h 2 )*, (5)
where m c is the effective mass of an electron in silicon and h is Planck's
constant. 6
A similar expression for the distribution of holes is
P(x) =( M Z) g(E)F h {E)dE, (6)
J — CO
where F h (E) = 1 - F,(E) ; substitution for g(E) and F h (e) yields
p(x) = N° f ™ <*■<*> ~ E)idE • (7)
m) ff 'L 1 + exp «E fh - E)/kT) (7)
E v (x) is the valence band edge and E fh is defined as the pseudo-Fermi
level 7 for the holes. N° is the constant given by
N° = 47r(2m„/h 2 )*, (8)
where m v is the effective hole mass in silicon. 8 E c (x) and E v (x) are func-
tions of <t>i(x) given by:
E e (x) = E°- qfrix), (9)
E v (x) = E% - q<t>i(x). (10)
E° and E° c are the valence and conduction band edges at x = « , and
q is the magnitude of an electronic charge. p(x) and n(x) are functions
of the yet undetermined constants E f and Efh. Later they will appear
only in the forms
v = (E°- E,)/kT,
7,' = (E° - E fh )/kT.
Since only difference terms appear, no energy reference need be estab-
lished for the model, y will be determined from the electron charge
neutrality condition at x ™ oo . q' is fixed by the total amount of free
positive charge Q+ = f hl °° p(x)dx in the device. In reality, determina-
1012 THE BELL SYSTEM TECHNICAL JOURNAL, JULY-AUGUST 1973
tion of r\' is difficult, so the alternative scheme of choosing v' and cal-
culating the resulting Q+ is employed. Substituting eqs. (2), (4), and (7)
into eq. (1) results in the system of differential equations which charac-
terize the device described above. In the appendix these equations are
reduced in a straightforward manner by using the Boltzman approxi-
mation; justification for its use is given. The resulting equations which
will be solved are (' denotes d/dy) :
flXy) = 0, ^ y ^ h/\, (11)
tf(y) = a - dm,/mJ*e"r+*M, hi/\ ^ y ^ h*/\, (12)
t!(y) - e* tiv) - 1 - (m v /m e )*e'"-+ tW , h*/\ £ y ^ », (13)
where
4, = q <j>/kT,
y = x/\
x = VkZVgWjy,
N t = N°(kT)*,
N c = N° c (kT)i,
t a = the bandgap energy in kT's,
v = (E° — E f )/kT is a constant dependent
on doping levels,
v ' = (E° c — E fh )/kT is a parameter depend-
ing on the pseudo-Fermi level,
p = V + v' — e a ,
a = N A /N D .
m, and m e are the effective masses of holes and electrons in silicon;
they are 1.08m and 0.59 m , respectively. 9 Equations (11), (12), and
(13) satisfy boundary conditions
*i(0) = V (14)
and, consistent with eqs. (9) and (10),
^ 3 (°°)=0. (15)
At y = hi/\ and y = hi/\ the continuity conditions are
*i(Ai/X) = Hhi/X) (16)
Hh2/\) = ^(VM (17)
BURIED-CHANNEL CCD 1013
II. NUMERICAL SOLUTION
The system of differential equations (11), (12), and (13) together
with boundary equations (14), (15), (16), and (17) are solved numeri-
cally by the method of finite differences.
For computational convenience the length of the device will be
truncated from the whole half-line to the segment [0, L]; the distance
L is chosen such that | ih( °° ) — ^i{L) | and | ^ 3 '( <*> ) — &a(L) \ are
suitably small. [0, L] is partitioned by a mesh of N points and the
solution at each point y, is denoted by \f/\ The mesh lengths (distance
between successive points) are 5 , 8 P , and S„ in the oxide, p-region, and
n-region, respectively, and are chosen such that Vni = hi/X and
y n, = h 2 /\. The boundaries are ^° = ^i(O) = V„ and ip N = f 3 (L)
« ^ 3 (oo) =0. The subscripts on \f/(y) can now be dropped since the
superscripts identify the solution point. The second derivative is ap-
proximated by the second difference
<W?/) _ ,,,,, O.,,- , ,HW«
dy'<
= (+ i+l - 2f< + t*- l )/S*
in each region; 5 is one of 8 OI 8„, or 8 P as appropriate. Using this ap-
proximation to discretize eqs. (11), (12), and (13) results in the matrix
equation
MO<]= IXy«,*9] (18)
which has rows generated by eq. (19); p(yt, $*) is defined in the three
regions by the right-hand side of eq. (19).
f0, i < N h
p+i - ty + ^•- 1 = U<r - (m,/m e )*e'*-* t )8 v , N x <i< N a , (19)
[( e W _ 1 _ (m,/w» e )*e"«-**)fi BJ i > N 2 .
Rows corresponding to the solutions at the boundary points t/at, and
yti, are obtained from eqs. (16) and (17) where the first derivative is
approximated by
=Ff(y) « (+W(y) - W(V±8)
+ fy(y ± 25) - 2+(y ± 35))/65, (20)
and 8 is the appropriate mesh size for each region (the positive or for-
ward derivative uses points to the left of y). Equation (20) results from
simultaneous solution of the Taylor expansions for \f/(y ± 8), \f/(y ±25),
and \J/(y ± 35) with terms 0(5 4 ) and higher dropped; the result is ac-
curate to 0(5 3 ) which is consistent with the second difference approxi-
mation to d 2 \J//dy 2 . In the oxide layer, a first-difference approximation
fix) = (*» - **-*)/«. (21)
1014 THE BELL SYSTEM TECHNICAL JOURNAL, JULY-AUGUST 1973
is adequate since the solution in this region is linear. The boundary-
equations are:
for i = Ni,
-(6* p ci/*.€»)^ w »- 1 + (11 + 6W«o€2)* Wl
- W" + fy Nl+i - 2^+ 3 = 0; (22)
for i = N 2 ,
(-1/*p)(2^ jv »- 8 - 9^'- 2 + 18^- J ) + ll**»(l/*» + 1/« P )
- (1/S n )(18^ s+1 - fy N > +2 + 2^ Ar »+ 3 ) = 0. (23)
The matrix A in eq. (18) is tridiagonal except for the Nith. and JV 2 th rows.
For each choice of the parameter p , the [$*] of eq. (18) are solved for
iteratively using a combination of successive under- and over-relaxa-
tion 10 (SOR) on the equation
[>•'] = M-»[>&**0J ( 24 )
The procedure described here differs slightly from the usual SOR in
that the transformed equation (24) is solved instead of eq. (18). For
the t'th row of eq. (24), the j -f- 1th estimate of ^* is given as
fti = # + "(-A*' - #), (25)
where w is a relaxation parameter with values < o> ^ 2. $' is a New-
ton's method solution of eq. (26), a transcendental equation in ^result-
ing from the ith row of eq. (24) with the remaining N — 2 variables ^*
held constant. The coefficients Out 1 in eq. (26) are elements of A' 1 .
f - E a* l p{y k , ftO - a^piyi, f ) - "£ a£ l p(Vk, *5) = 0. (26)
A-l k-»+l
Criteria for convergence of the process, as well as the choice of the
value of w, is based on the residual rj defined at the jth. iteration as
r,= Z |*}-fti I • (27)
It can be shown that if (18) were a linear system of equations and if f
were the true solution at ?/, then
sup If -ft.il/aq> \?\*CMr M (28)
where C(co) is a constant dependent on the choice of »." For the optir
mum value of v, C(oj) is 0(N) while for values of cu only slightly different
BURIED-CHANNEL CCD 1015
from the optimum C(oj) can be 0(iV 2 ). In the computations presented
in the next section, N ^ 100 so the iterative process was stopped when
the residual terms were less than 10 -8 ; this allowed margin for the fact
that (18) is nonlinear. A discussion of the choice of to is necessarily even
more heuristic. It was found that for a mesh of N = 50 points and large
negative values of p corresponding to small positive charge densities
(i.e., e"»~*'' « 1 for all i), the iterative scheme was convergent for any
choice of cu and any initial estimate of the {^*} . For less-negative values
of p and the corresponding larger values of e Po "~*', the scheme was con-
vergent for over-relaxation (w > 1) only if r, < 10 for all j (approxi-
mate figure) ; but using under-relaxation the process was well behaved
with r/s as great as 10 4 . It should be pointed out that due to the ex-
ponential nature of the right-hand side of eq. (19) even a very good initial
estimate typically resulted in ri > 10 3 when over-relaxation was ap-
plied. Although the SOR process described above may well be stable 12
regardless of r, values, the magnitude of the exponential terms in (19)
limits machine computations. The combination of under- and over-
relaxation detailed below eliminates this problem.
Good initial estimates for the {^'} were obtained by choosing the p„
Values with equal spacing Ap; the first p value being the small positive
charge-density case described above. With Ap < 10 a linear estimate
of the {$'] Tor successive values of p was adequate. The iterations were
under-relaxed (« = 0.5) until the residual term was less than 10; suc-
cessive iterations were over-relaxed so as to increase the rate of con-
vergence. The choice of as for the SOR steps was not critical; values in
the range 1 < w 2g 1.5 had approximately the same rate of conver-
gence. Values above 1.5 did not converge for all values of p .
Total positive charge in the device for a given value of p is calcu-
lated by
Q + = f ""' e'°-**Mdy + r 6 e"-**^dy. (29)
Using a piecewise linear approximation for \J/(y) on [ijk, Vk+\],
+(y) SS P + ((» - ¥*)/*)(#*» - **), (30)
and summing the integrals over each such interval in [_Vn v L~] gives the
approximation
Q + ^S P £ e-ie-*' - e-* i+l )/(f i+1 - *»')
+ S n Z e*>(e-*< - e-*' + W' +1 - M, (31)
1016 THE BELL SYSTEM TECHNICAL JOURNAL, JULY-AUGUST 1973
which can be modified to include the case \f/ i+1 = ip\ Dimensional
charge in coulombs can be calculated from Q+ by Q = Q+Nd^-
The operational scheme may be summarized as follows:
(i) Choose a value of p corresponding to small total charge Q +
( Po -» - * in eq. (24) => Q+ -^ 0),
(it) Solve for \f/(y), thus determining p(y),
(Hi) Integrate p(y) to find Q+,
(iv) Increment p by Ap and repeat (ii), (in), and (iv).
The technique of starting with Q+ « and slowly adding positive
charge to the device is important since it is this scheme of operation,
in conjunction with the successive under- and over-relaxation, that
avoids the exponential overflow limitations in machine computation.
III. COMPUTATIONAL RESULTS
In this section, numerical results for two specific device configura-
tions will be given. The first device has a constant doping profile
(Na = a constant) and dimensions
hi/X = 0.48,
h 2 /\ = 12.48,
L = 42.48.
4 8 12 16 20 24 28 32 36 40 44
x/\- DISTANCE IN DEBYE LENGTHS
Fig. 2 — Potential versus distance for a buried-channel device with uniformly
doped p-region and V„ = — 4 volts.
BURIED-CHANNEL CCD
1017
18
16
(9 6
2 -
4 8 12 16 20 24 28
x/a- DISTANCE IN DEBYE LENGTHS
Fig. 3 — Charge density versus distance for a buried-channel device with uniformly
doped p-region.
X in this case is approximately 0.415 /mi. Doping levels are
N A = 2 X 10 15 cm" 3 ,
JV D = 1X 10 14 cm" 3 ,
so
a = Na/Nd = 20.
N e is calculated by 13
N c = 4.831 X 10"(m«/m.)»r»,
with T = 300° K and m c as before.
Using eq. (46), 77 is found to be
77 = 12.09.
The nondegeneracy condition for this device may now be stated using
eq. (45)
6, - 17' + +(y) > 3.5,
1018 THE BELL SYSTEM TECHNICAL JOURNAL, JULY-AUGUST 1973
-200
-250 -
u. (/I
O
-300 -
z £ -350 -
-400 -;
20
40
60 80 100
TOTAL CHARGE
Fig. 4 — Potential at the Si-Si0 2 interface versus total positive charge in a buried-
channel device with a uniform p-region and V = — 4 volts.
or
V' ~ e g - +(V) < -3-5.
Adding rj to both sides of the inequality gives
Po - Hy) < V - 3.5.
-800
12 16 20 24
x /A -DISTANCE IN DEBYE LENGTHS
28
32
36
Fig. 5 — Potential versus distance for a buried-channel device with a Gaussian
p-region doping profile and V = —4 volts.
BURIED-CHANNEL CCD
1019
36 -
o
o 3? ~
o
u
24 -
x 2C -
t
2 16
g 12
<
8 -
- K-\Q t = 187.6
| *V— 135.0
*-V -965
-
K-4- 64.4
-
(*-*- 36.6
AilW
Oj
I
1 .
4 8 12 16
x/X- DISTANCE IN DEBYE LENGTHS
20
Fig. 6 — Charge density versus distance for a buried-channel device with a Gaussian
doping profile.
Then for each choice of the parameter p , the solution set {^'} must
satisfy:
Po - +< < 8.59. (32)
From eq. (12) it is clear that if one assumes the Boltzmann approxi-
mation to hold, then as long as Q + < <r(h 2 — hi)/X the equilibrium con-
dition causes the lineal charge density to have constant sign so
Po — f* < log <j(m c /m v )*
(33)
for all i; for this device, eqs. (32) and (33) are always consistent if
N A < 1.3 X 10 18 cm- 3 .
Computation was performed using a 90-point mesh and took ap-
proximately 4.8 minutes of processor time on a Honeywell 6000-series
machine.
Figures 2, 3, and 4 summarize the computational results. Figure 2
shows potential solutions {if/*} versus the points ?/, for a family of p„
1020 THE BELL SYSTEM TECHNICAL JOURNAL, JULY-AUGUST 1973
£
-200
-i
c
>
J
<
5
-250
—
EC
UJ
I
H
Z
_1
-300
<
1-
2
UJ
K
O
Q-
-350
UJ
y
<
u.
BE
"max = 46
LU
1-
V Q = -4
z
-400
— ^"^
o
in
Ch
-450
1 1
1 1 1 1 1
40 60 80 100 120 140 160
TOTAL CHARGE
Fig. 7— Potential at the Si-Si0 2 interface versus total positive charge in a buried-
channel device with a Gaussian p-region doping profile and V a = — 4 volts.
values with V = — 4 volts. The range of total charge values is indi-
cated. Figure 3 is a plot of the linear charge density with Q+ values
indicated. Figure 4 is Q+ plotted against $ Nl (the oxide-interface
potential).
A realistic modification to the device studied thus far is to allow a to
have y variation. The second set of results presented are for the same
device as described above but with the p-region having a doping profile
a(y) = (2>. + l)tf-l<r-»i>/(ta-ii)l« ^ (*>■+» _ 2. (34)
The values for h h h 2 , and V are as before, L = 32.48, and D, = 46.
[This value of D, corresponds to an average doping level in the p-layer
of N A = 1.9 X 10 16 /cm 8 , and eq. (34) describes a doping profile as if
the p-layer were formed by drive-in diffusion. 14 ] The solutions \\l/ 1 }
must still satisfy eq. (32). Figures 5, 6, and 7 are a summary of results.
Comparison of the two cases shows that for the doping profile in eq.
(34) the potential minimum is greater and the "channel" is shifted
toward the oxide. In both cases the added positive charge is contained
entirely in the p-region and it fills the region starting from the side
remote from the oxide interface (see Figs. 2 and 5).
IV. ACKNOWLEDGMENTS
The work presented here has benefitted from constructive sugges-
tions by N. L. Schryer, G. E. Smith, R. J. Strain, and D. A. Kleinman.
BURIED-CHANNEL CCD 1021
Special thanks goes to J. McKenna for his many contributions during
the course of this work and for his comments regarding the preparation
of this paper.
APPENDIX
Substituting eqs. (2), (4), and (7) into eq. (1) results in (' denotes
d/dx) :
4>'i(x) = 0, 0gx£h h (35)
Mx)--]N A -N.J_ m 1 + exp L{Efh _ E)/kT1 j- ,
h ^ x ^ hi, (36)
»/ x g|»„r (E-E e (x))*dE
m f '•« (E 9 (x) - E)»dE \
"•}-* 1 + exp [(£,„-£) AT] J'
ht ^ x S w . (37)
Transform the integral involving E c (x) by making the substitutions:
e = (E - E c (x))/kT,
E c (x) = E° c - q<p{x),
v = (El - E f )/kT.
For the integrals involving E v (x) make the substitutions:
e = (E v (x) - E)/kT,
E,(x) = E% - q V (x),
rj' = (E° c - E fh )/kT,
e = (E° - E°)/kT.
For both cases let
4>{x) = q<p(x)/kT
and make the change of variables
V = */\
where
x = VkTVgWD.
These substitutions and some straightforward manipulation reduce
eqs. (35), (36), and (37) to
tib) =0, ^ y g hi/\, (38)
,»(»\ __N A _N% n-n,/"* e^de
* 2{y) - N D N D {kI) in 1 + exp (e + «. - r,' + fcfo)) '
fej/X ^ y ^ A 2 /X, (39)
1022 THE BELL SYSTEM TECHNICAL JOURNAL, JULY-AUGUST 1973
tfdr) = jt d (kT)*J o 1 + exp p + - _ ^ |(y) -, - 1
tf D ^ J 7o 1 + exp (£ + «.- f,' + **(»)) '
A2A ^ 1/ ^ ». (40)
These equations are further simplified by making the classical Boltz-
mann approximation for the electrons :
r *** ~r ^e (41)
Jo 1 + exp (e + v ~ t(.y)) Jo exp(e + 77.- $(y))
which holds since y - \p(y) is always positive by several kT (see Fig.
8). The last result is integrable and simplifies the approximation to
T(f) exp (+(y) ~ v), (42)
where r(-) is the usual gamma function. This approximation has been
shown to be less than one percent in error if 16
t, - +(y) > 3.5. (43)
Similarly, for the holes,
/
€^€
^ r(|) exp w - i„ - +(y)) (44)
o 1 + exp (e + e - y' + *(y))
if
e B ~ v' + *(V) > 3.5. (45)
The validity of this last approximation is not clear since +(y) is de-
pendent on i}' and \p(y) becomes large and negative. Note that requiring
the inequalities (43) and (45) to hold is equivalent to requiring the
device to be nondegenerate. Computational results show this approxi-
mation to be consistent.
77 is expressed in terms of the constants N c , N v , and N D by requiring
electron charge neutrality at x = -» . It should be stated that the cor-
rect condition to use at this point is complete charge neutrality. This
requirement, however, is computationally indistinguishable from the
algebraically simpler condition used here. Using eqs. (40) and (15), the
neutrality condition is
(N e /N D )e-"T(%) -1=0,
so
N C /N D = e'/m). (46)
BURIED-CHANNEL CCD
1023
Fig. 8 — Band diagram for buried-channel device.
Equation (46) defines the Fermi level for the electrons. Using eqs. (5)
and (8),
N V /N D = (N,/N C )-(N C /N D ) = (m./w.)Wr(t), (47)
where m v and m c are the effective masses of holes and electrons in
silicon; they are 1.08 m and 0.59 m , respectively. 10 Using eqs. (42),
(44), (46), and (47) in eqs. (38), (39), and (40) results in equations
(11), (12), and (13) which are to be solved.
REFERENCES
1. Walden, R. H., Krambeck, R. H., Strain, R. S., McKenna, J., Schryer, N. L.,
and Smith, G. E., "The Buried Channel Charge Coupled Device," B.S.T.J.,
61, No. 7 (September 1972), pp. 1635-1640.
2. Amelio, G. F., "Computer Modeling of Charge-Coupled Device Characteristics,"
B.S.T.J., 51, No. 3 (March 1972), pp. 705-730.
3. McKenna, J., and Schrver, N. L., "The Potential in a Charge Coupled Device
With No Mobile Mmority Carriers and Zero Plate Separation," B.S.T.J.,
52, No. 5 (May-June 1973), pp. 669-696.
4. McKenna, J., and Schryer, N. L., "The Potential in a Charge-Coupled Device
With No Mobile Minority Carriers," unpublished work.
5. Strain, R. J., and Schryer, N. L., "A Nonlinear Diffusion Analysis of Charge-
Coupled-Device Transfer," B.S.T.J., 50, No. 6 (July- August) , pp. 1721-1740.
6. Blakemore, J. S., Semiconductor Statistics, New York: Pergamon Press, 1962,
pp. 75-77.
1024 THE BELL SYSTEM TECHNICAL JOURNAL, JULY-AUGUST 1973
7. Grove, A. S., Physics and Technology of Semiconductor Devices, New York : John
Wiley and Sons, 1967, p. 162.
8. Blakemore, J. S., op. cit., p. 81.
9. Blakemore, J. S., op. cit., pp. 58-63. , XT .
10. Varga, R. S., Mafria: Iterative Analysis, Englewood Cliffs, N. J. : Prentice-Hall,
1962, p. 59.
11. McKenna, J., unpublished work. ,
12. Richtmeyer, R. D., and Morton, K. W., Difference Methods for Initial Value
Problems, 2nd Ed., New York: Interscience, 1970, p. 9.
13. Blakemore, J. S., op. cit., p. 82.
14. Grove, A. S., op. cit., p. 43.
15. Blakemore, J. S., op. cit., p. 358.